###########################################
#REPLICATION FILE FOR VOX POPULI, VOX DEI?
#by ADAM RAMEY
#SECTION 6 (incl. FIGURES 5, 8, 9, and 10)
###########################################

library(basicspace)
library(foreign)
library(ggplot2)
library(doParallel)
library(scales)
library(plyr)
library(MASS)
library(ordinal)
library(reshape2)
library(texreg)

##be sure to change the line below so as to indicate where the data and scripts are
setwd("~/Downloads/Replication Materials/")

#this commented out line has the results all done; uncomment if you want to skip the
#estimation
#load("GovData.Rda")

ds <- read.dta("cces_cumul.dta")
ds <- subset(ds, ds$year %in% seq(2009,2012,by=1))

attach(ds)

##national actors
self <- ideo7
obama <- ideo7_obama
democrats <- ideo7_dems
republicans <- ideo7_gop
scotus <- ideo7_scotus
teaparty <- ideo7_teaparty

placement_matrix <- cbind(self, obama, democrats, republicans)

#take care of missing values
for (j in 1:ncol(placement_matrix)){
  placement_matrix[,j] <- as.numeric(placement_matrix[,j])
  placement_matrix[,j] <- replace(placement_matrix[,j], placement_matrix[,j]==8 | is.na(placement_matrix[,j]), 999)
}  


overall <- aldmck(placement_matrix, polarity=2, respondent=1, missing=999, verbose=TRUE)

#make a quick density plot (FIGURE 5)
g1 <- ggplot(overall$respondents, aes(idealpt)) + 
  geom_density(fill="grey",colour="grey") +
  xlab("Respondent Ideal Points") + 
  ylab("Density") +
  theme_bw()

ggsave(g1, file="voterdensity.pdf",
       height=8, width=6)

##Now, we estimate governors
states <- names(tapply(state_post,state_post,length))[!is.na(tapply(state_post,state_post,length))]


results <- NULL

for (k in 1:length(states)){
  
  for (yr in seq(2009,2012,by=1)){
    
    subsetting <- (state_post %in% states[k] | state_pre %in% states[k]) & year %in% yr
    tmp <- placement_matrix[subsetting,]
    governor <- replace(cbind(ideo7_gov), 
                        ideo7_gov==8 | is.na(ideo7_gov), 999)[subsetting]
    pm <- cbind(tmp, governor)
    scaling <- aldmck(pm, polarity=2, respondent=1, missing=999, verbose=FALSE)
    
    registerDoParallel(cores=8)
    scaling_boot <- foreach(i=1:200, .combine='cbind') %dopar% 
      aldmck(pm[sample(1:nrow(pm),nrow(pm),replace=TRUE),], polarity=2, respondent=1, 
             missing=999)$stim
    
    slope <- (overall$stim[3]-overall$stim[2])/(scaling$stim[3]-scaling$stim[2])
    rescaling <- slope*scaling$stim - slope*scaling$stim[2] + overall$stim[2]
    
    rescaling_boot <- slope*scaling_boot[4,] - slope*scaling$stim[2] + overall$stim[2]
    
    
    #save results for elites
    tmp2 <- c(rescaling[4], quantile(rescaling_boot,c(0.05,0.95)))
    tmp2 <- data.frame(names(rescaling)[4],t(tmp2))
    colnames(tmp2) <- c("governor","idealpt","lower","upper")
    
    #voters
    ip <- overall$resp$idealpt
    pid <- as.numeric(ds$pid7)
    
    voter_boot <- foreach(i=1:200, .combine='cbind') %dopar% 
      mean(ip[subsetting][sample(1:nrow(pm),nrow(pm),replace=TRUE)],na.rm=TRUE)
    
    rep_boot <- foreach(i=1:200, .combine='cbind') %dopar% 
      mean(ip[subsetting & pid %in% c(6,7)][sample(1:nrow(pm),nrow(pm),replace=TRUE)],na.rm=TRUE)

    dem_boot <- foreach(i=1:200, .combine='cbind') %dopar% 
      mean(ip[subsetting & pid %in% c(1,2)][sample(1:nrow(pm),nrow(pm),replace=TRUE)],na.rm=TRUE)
    
    #put together    
    res_tmp <- data.frame(states[k], yr, 
                          mean(ip[subsetting],na.rm=TRUE),
                          quantile(voter_boot,c(0.05),na.rm=TRUE),
                          quantile(voter_boot,c(0.95),na.rm=TRUE),
                          mean(ip[subsetting & ds$pid3=="Republican"],na.rm=TRUE),
                          quantile(rep_boot,c(0.05),na.rm=TRUE),
                          quantile(rep_boot,c(0.95),na.rm=TRUE),
                          mean(ip[subsetting & ds$pid3=="Democrat"],na.rm=TRUE),
                          quantile(dem_boot,c(0.05),na.rm=TRUE),
                          quantile(dem_boot,c(0.95),na.rm=TRUE),
                          tmp2[-1])
    
    colnames(res_tmp) <- c("state","year","voters","voters_low","voters_high",
                           "rep","rep_low","rep_high",
                           "dem","dem_low","dem_high",
                           "governor","gov_low","gov_high")
    
    results <- rbind(results, res_tmp)
    
    cat(colnames(res_tmp), "\n", as.matrix(res_tmp), "\n\n")
  }
}

results <- results[results$state!="District of Columbia",]

results$rep[is.nan(results$rep)] <- ((results$rep_low+results$rep_high)/2)[is.nan(results$rep)]
results$dem[is.nan(results$dem)] <- ((results$dem_low+results$dem_high)/2)[is.nan(results$dem)]

##merge the results
ds$voter_ip <- overall$respondents$idealpt
ds$stateyear <- paste(ds$state_post,ds$year)
results$stateyear <- paste(results$state,results$year)
dd <- merge(ds, results, by.x="stateyear", by.y="stateyear")


##recode the governor approval variable
dd$app_gov <- as.numeric(dd$approval_gov)
dd$app_gov <- (dd$app_gov==1)*2 + (dd$app_gov==2)*1 + (dd$app_gov==5)*0 +
  (dd$app_gov==3)*-1 + (dd$app_gov==4)*-2

dd$gov_voter_dist <- abs(dd$voter_ip-dd$governor)
dd$sameparty <- ((dd$pid3=="Republican" & dd$governor_party=="Republican")|
                 (dd$pid3=="Democrat" & dd$governor_party=="Democrat"))*1   

mod <- (glm(app_gov>0 ~ gov_voter_dist + sameparty, data=dd, family=binomial))

#make predicted values
preds <- predict(mod, 
                 newdata=data.frame(gov_voter_dist=rep(seq(0,2,by=0.05),2),
                                    sameparty=c(rep(1,length(seq(0,2,by=0.05))),
                                                rep(0,length(seq(0,2,by=0.05))))), 
                 type="response",
                 se.fit=TRUE)

#make the ggplot by melting predictions and errors
newdata <- data.frame(gov_voter_dist=rep(seq(0,2,by=0.05),2), 
                      sameparty=c(rep(1,length(seq(0,2,by=0.05))),
                                  rep(0,length(seq(0,2,by=0.05)))),
                      se=preds$se.fit, 
                      fits=preds$fit)

#newdata <- melt(newdata, id.vars = c("se","fits","sameparty"))

##Figure 10
g2 <- ggplot(newdata, aes(gov_voter_dist,fits,fill=factor(sameparty),linetype=factor(sameparty))) + 
  geom_line() + 
  geom_ribbon(alpha=0.3, aes(xmin=gov_voter_dist,xmax=gov_voter_dist,
                             ymin=fits-1.96*se,ymax=fits+1.96*se, fill=factor(sameparty),
                             linetype=NA)) +
  scale_fill_grey(name = "", 
                  breaks=c("1","0"), labels=c("Same Party", "Different Party"),
                  start = 0.0, end = 0.6) +
  scale_linetype(name = "", 
                 breaks=c("1","0"), labels=c("Same Party", "Different Party")) +
  xlab("Distance between Voter and Governor Ideal Point") + 
  ylab("Pr(Strongly or Somewhat Approve of Governor)") +
  theme_bw()

ggsave(g2, file="GovApp.pdf",
       height=6,width=9)


###make a graph of the median voter over time
res2 <- data.frame(rep(results$state,4), 
                   rep(results$year,4), 
                   c(results$governor,results$voters,results$rep,results$dem),
                   c(results$gov_low,results$voters_low,results$rep_low,results$dem_low),
                   c(results$gov_high,results$voters_high,results$rep_high,results$dem_high),
                   rep(c("Governor","Mean Voter","Republican Voters","Democratic Voters"), each=nrow(results)))

colnames(res2) <- c("state","year","idealpts","lower","upper","who")

##add variable for switching
switch_data <- subset(res2,res2$who=="Governor")
switch_data$change <- 0

for (i in 2:nrow(switch_data)){
  if (switch_data$state[i]==switch_data$state[i-1] & 
      abs(switch_data$ideal[i]-switch_data$ideal[i-1])>0.3) switch_data$change[i] <- 1
}

switch_data <- switch_data[switch_data$change==1,]
switch_data$switch_year <- switch_data$year

res2 <- merge(switch_data,res2,all=TRUE)

##Figure 8
g3 <- ggplot(res2,aes(year,idealpts,group=who,fill=who,linetype=who)) + 
  geom_ribbon(aes(xmin=year,xmax=year,ymin=lower,ymax=upper,linetype=NA),alpha=0.3) +
  geom_line(aes(year,idealpts)) +
  geom_point(aes(year,idealpts,shape=who)) +
  geom_vline(aes(xintercept=switch_year), linetype="dotted") +
  facet_wrap(~state,ncol=10) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  scale_fill_grey(name = "", 
                  labels=c("Democratic Voters", "Governor", "Mean Voter", "Republican Voters"),
                  start = 0.0, end = 0.7) +
  scale_colour_grey(name = "", 
                    labels=c("Democratic Voters", "Governor", "Mean Voter", "Republican Voters"),
                    start = 0.0, end = 0.7) +
  scale_linetype(name = "", 
                 labels=c("Democratic Voters", "Governor", "Mean Voter", "Republican Voters")) +
  scale_shape(name = "", 
              labels=c("Democratic Voters", "Governor", "Mean Voter", "Republican Voters")) +
  #scale_fill_manual("", values = alpha(c("blue", "black", "purple", "red"), 0.3)) +
  #scale_colour_manual("", values = alpha(c("blue", "black", "purple", "red"), 0.9)) +
  xlab("Year") + ylab("Ideal Points") 

ggsave(g3, file="govs.png", height = 10, width = 14) 

##Figure 9
g4 <- ggplot(results,aes(year,abs(voters-governor))) + 
  geom_line() +
  facet_wrap(~state,ncol=10) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + 
  xlab("Year") + 
  ylab("Ideological Divergence between Governors the Voters") 

ggsave(g4, file="govs_div.png", height = 10, width = 14) 



##look at the Peress-Batista-Richman stuff
mod1 <- lm(governor ~ voters, data=results)
mod2 <- lm(governor ~ voters, data=subset(results,year==2009))
mod3 <- lm(governor ~ voters, data=subset(results,year==2010))
mod4 <- lm(governor ~ voters, data=subset(results,year==2011))
mod5 <- lm(governor ~ voters, data=subset(results,year==2012))

texreg(list(mod1,mod2,mod3,mod4,mod5))



